Analysis of single nuclear chromatin accessibility reveals unique myeloid populations in human pancreatic ductal adenocarcinoma

Abstract Background A better understanding of the pancreatic ductal adenocarcinoma (PDAC) immune microenvironment is critical to developing new treatments and improving outcomes. Myeloid cells are of particular importance for PDAC progression; however, the presence of heterogenous subsets with different ontogeny and impact, along with some fluidity between them, (infiltrating monocytes vs. tissue‐resident macrophages; M1 vs. M2) makes characterisation of myeloid populations challenging. Recent advances in single cell sequencing technology provide tools for characterisation of immune cell infiltrates, and open chromatin provides source and function data for myeloid cells to assist in more comprehensive characterisation. Thus, we explore single nuclear assay for transposase accessible chromatin (ATAC) sequencing (snATAC‐Seq), a method to analyse open gene promoters and transcription factor binding, as an important means for discerning the myeloid composition in human PDAC tumours. Methods Frozen pancreatic tissues (benign or PDAC) were prepared for snATAC‐Seq using 10× Chromium technology. Signac was used for preliminary analysis, clustering and differentially accessible chromatin region identification. The genes annotated in promoter regions were used for Gene Ontology (GO) enrichment and cell type annotation. Gene signatures were used for survival analysis with The Cancer Genome Atlas (TCGA)‐pancreatic adenocarcinoma (PAAD) dataset. Results Myeloid cell transcription factor activities were higher in tumour than benign pancreatic samples, enabling us to further stratify tumour myeloid populations. Subcluster analysis revealed eight distinct myeloid populations. GO enrichment demonstrated unique functions for myeloid populations, including interleukin‐1b signalling (recruited monocytes) and intracellular protein transport (dendritic cells). The identified gene signature for dendritic cells influenced survival (hazard ratio = .63, p = .03) in the TCGA‐PAAD dataset, which was unique to PDAC. Conclusions These data suggest snATAC‐Seq as a method for analysis of frozen human pancreatic tissues to distinguish myeloid populations. An improved understanding of myeloid cell heterogeneity and function is important for developing new treatment targets in PDAC.

1. snATAC-Seq can be used to assay frozen benign pancreatic and PDAC tissues.
2. snATAC-Seq data can be used to determine the regulome and molecular characteristics of cell populations.

Myeloid cells from the PDAC TME can be characterized using snATAC-Seq. INTRODUCTION
Although significantly improved over the past decade, the 5-year survival rate for pancreatic ductal adenocarcinoma (PDAC) is only 12%.The unique stromal makeup of the PDAC tumour, including activated fibroblasts and suppressive immune cells, is a critical component driving tumour progression and resistance to therapy.To enhance therapeutic options for PDAC patients, a better understanding of this unique tumour microenvironment (TME) is needed.
Immune cell infiltration of the PDAC TME is an important subject of inquiry as immunotherapies with success in other cancer types have failed in PDAC. 1 In particular, the myeloid component of the PDAC TME has been identified as having a critical role in this immunosuppression and treatment failure.Circulating monocytes were correlated with a worse overall survival in PDAC patients 2 ; however, clinical trials blocking monocyte recruitment into the tumour combined with FOLFIRINOX 3 or gemcitabine/nab-paclitaxel 4 demonstrated limited efficacy.Another population of myeloid cells, tissue-resident macrophages (TRMs), are present before tumour initiation and are the prominent source of immunosuppressive, pro-fibrotic macrophages in the PDAC TME in murine models. 5Meanwhile, M1 and M2 remain the most common characterisation of macrophage populations in the TME, but individual macrophages exist along more of a continuum than this nomenclature reflects. 6Such studies demonstrate a need for refined characterisation of myeloid source and function in the PDAC TME to understand pathogenesis and uncover novel therapeutic strategies.
Recently, single cell sequencing technologies have been harnessed to determine cell function and composition in the TME.][9][10][11][12] Since stroma can comprise much of the tumour volume, 8,13 scRNA-Seq has been employed to identify the immune cell milieu and ascertain how these cells influence tumour progression, 10 patient prognosis 7,11,12 and chemotherapy response. 11Given the importance of myeloid cells in PDAC progression, 2,5 scRNA-Seq studies have also explored the phenotype and function 11,12 as well as source, 12 of macrophages.Although the transcriptome of cells within the PDAC TME has been heavily studied, 11,[14][15][16] there has been limited progress on the regulome of cells within the TME.To answer such questions, a method for analysing open chromatin, termed assay for transposase accessible chromatin sequencing (ATAC-Seq), can explore the epigenetic underpinnings of diseases. 17,18Open chromatin can also be utilised for cell identification due to its role in regulation of gene expression. 18,19A limited number of studies exploring single cell chromatin accessibility from freshly isolated PDAC specimens have shown an important role for KRAS regulation of enhancer activity 20 and dysregulation of methylation 7 in PDAC progression, demonstrating the importance of epigenomics in understanding pancreatic diseases.Furthermore, open chromatin and transcription factor activity can be used to determine myeloid cell origin, 19 which is relevant for understanding the highly plastic myeloid populations that are present in the PDAC tumour and the role for these cells in cancer progression.
In the present study, we applied single nuclear ATAC-Seq (snATAC-Seq) to frozen benign and PDAC patient pancreatic specimens.Our results justified snATAC-Seq as an effective method for retrospective analysis of tumours to identify unique cell populations within the TME.Through snATAC-Seq, we combined information on transcription factor binding and gene promoter activity to identify cell populations and infer their functions in benign pancreatic and PDAC samples.Utilising these findings, we then explored the impact of distinct myeloid cell populations on patient survival.

Patient tissue collection
Prior to sample collection, Institutional Review Board approval was obtained from West Virginia University (#1903496995) and all patients signed informed consent.
All patients underwent surgical resection with pancreaticoduodenectomy as part of standard of care.Benign and PDAC samples for snATAC-Seq were treatment naïve.Tissue was obtained from the pathology frozen section suite immediately upon surgical resection.Diagnosis was confirmed by pathology.PDAC tumour tissue or benign pancreatic tissue was frozen alone (nucleus isolation) or in OCT (optimal cutting temperature compound; immunofluorescence staining) and stored at −80 • C.

Immunofluorescent staining
Sections (10 μm) of OCT-embedded frozen tissues were fixed in 10% formalin then permeabilised with .1% Triton-X for 10 min and blocked for

Nuclei isolation from frozen tissue
Isolation procedure was adapted from nuclei isolation from complex tissues for Single Cell Multiome ATAC + Gene Expression Sequencing (CG000375).Nuclei were isolated from 50 mg of frozen tissue cut into five small sections with a razor blade.The tissue pieces were homogenised for 1.5 min with a tissue homogeniser and pestle in 10 mM Tris-HCl, 10 mM NaCl, 3 mM MgCl 2 , .1% Nonidet P40 substitute and 1 mM dithiothreitol (DTT) in nuclease-free water.Tissue was dissociated in lysis buffer for 5 min before passing through a 70 μm cell strainer.The solution was then centrifuged (500 × g, 4 • C, 5 min).The supernatant was removed and replaced with 1 mL 1% BSA in DPBS and allowed to incubate for 5 min before resuspension.The nuclei were centrifuged again (500 × g, 4 • C, 5 min).The pellet was resuspended in 1 mL 1% BSA in DPBS and passed through a 30 μm cell strainer.All steps were performed on ice.Nuclei were stained with 7-aminoactinomycin D (7-AAD) readymade solution 1:100 (Sigma, SML-1633) and sorted with a FACSAria III Cell Sorter.Nuclei were collected into 10% BSA.For snMultiomics, 1 U/mL RNase inhibitor (Sigma-Aldrich, 03335402001) was added to all final solutions.

Bulk ATAC-Seq
The Omni-ATAC Protocol 18 was followed.Libraries were prepared with sorted nuclei with a Tagment DNA Enzyme and Buffer Large Kit (Illumina, 20034198) and run on 2% agarose gel to observe the nucleosome phasing pattern and size select.The gel was cut for an approximately ∼200−600 bp band and subjected to gel elution using the Qiagen gel MinElute Gel Extraction Kit (Qiagen, 28606).

snATAC-Seq data pre-processing and clustering analysis
snATAC-Seq reads were aligned to the GRCH38(hg38) reference genome.Fragment counts were quantified using cellranger-atac count (version 2.1.0).Seurat 21 (version 4.2.1) and Signac 22 (version 1.7.0) were used for downstream analysis.Quality control was based on library features with settings: 1000 < nCount_ATAC (number of fragments) < 20 000, 500 < nFeature_ATAC (number of peaks) < 10 000, blacklist_fraction < 5%, nucleosome_signal < 2, TSS.enrichment > 2 and pct_reads_in_peaks (fraction of fragments within ATAC-Seq peaks) > 15%.Nucleosome_signal (strength of the nucleosome signal), TSS.enrichment (transcription start site [TSS] enrichment score as defined by ENCODE) and blacklist_fraction (ratio reads in genomic blacklist regions as defined by ENCODE) were calculated by the function NucleosomeSignal, TSSEnrichment and Frac-tionCountsInRegion of Signac, respectively.RunTFIDF was used to compute term-frequency inverse-documentfrequency normalisation on the matrix.FindTopFeatures was used to find the most frequently observed features.RunSVD was used to find the top 50 largest singular value and corresponding singular vectors of the matrix.RunHarmony 23 was used to correct batch effects.RunUMAP was used for dimensional reduction and visualisation via uniform manifold approximation and projection (UMAP).FindNeighbours was used to compute the nearest neighbours for the object.FindClusters was used to identify clusters of cells by a shared nearest neighbours modularity optimisation based on Smart Local Moving (SLM) clustering algorithms 24 with a resolution of .2. FindAllMarkers was used to identify the differentially accessible regions for each cluster by against all other clusters, by setting 'LR' as the test used, peak_region_fragments as the test variable, and logfc.threshold(fold change in log2 scale) as .25.
The RNA part of the snMultiome was used to annotate cell type and the ATAC part for validation of the cell type annotation in the snATAC-Seq (next section).For the RNA part, Seurat 21 (version 4.2.1) was used for downstream analysis.NormalizeData was used for normalisation and set normalisation method as 'CLR' (centered log ratio transformation).Cell cycle phase score was calculated by CellCycleScoring.Variable features were identified by FindVariableFeatures and the top 2000 were selected.Data were scaled and centered by ScaleData function and regress out the cell cycle's effect.RunPCA was used for principal component analysis (PCA) dimensionality reduction.RunHarmony 23 was used to integrate the dataset.RunUMAP and FindNeighbours were called as in the previous section.FindClusters was used to identify clusters of cells by a shared nearest neighbours modularity optimisation based on original Louvain clustering algorithms with a resolution of .2. Cell clusters annotated with marker genes are acinar cell (CUZD1 25 ), acinar to ductal metaplasia (ADM) cell (ONECUT1 26 ), ductal cell (KRT19 11,15 ), endocrine cell (PAX6 27 ), endothelial cell (TEK 28 ), fibroblast (COL1A2 29 ), myeloid cell (CD163, 12 CSF2RA 30 ) and T cell (IL7R 11 ) (Table S1).After annotation, an extended list of marker genes was predicted for each cell type against all other cell types through FindAllMarkers with 'wilcox' test and logfc.thresholdas .25 (Table S1).

Cell type annotation and validation for snATAC-Seq
Cell type annotation for the snATAC-Seq datasets was done through an integrative analysis with annotated cells from the snMultiome dataset.After differential accessible regions were identified for each cluster from the snATAC-Seq dataset (Table S2), promoters overlapping with these peaks were used to define signature genes for each cluster.Signature genes for each cluster from the snATAC-Seq dataset were then compared to marker genes for each cell type defined from the RNA part of the snMultiome data set for an overrepresentation analysis (observed number/expected number) (Figure S5D).For each snATAC-seq cluster, the cell type with the highest overrepresentation was assigned to that cluster.
We validated our snATAC-Seq annotations by four methods (see Section 3).First, expression of representative genes in snATAC-Seq clusters was examined in the snMultiome dataset to visually validate cell type annotation.Second, we identified the overrepresented molecular pathways for snATAC-Seq clusters with Gene Ontology (GO) enrichment of signature genes.Third, we aligned the snATAC-Seq data with the ATAC part of our snMultiome.Cell type annotation for the ATAC part of the snMultiome data was directly transferred from the RNA part of the snMultiome data.We applied the RunHarmony function in the Harmony package 23 for this alignment, visualised data in UMAP, projected their respective cell type annotation and visually examined their positional relationship.
Lastly, four scATAC-Seq datasets 20 were downloaded from Gene Expression Omnibus (GSE147726) (one patient tumour and adjacent nonmalignant tissue, another tumour and a histologically normal pancreas).While cell types are listed by Li et al., 20 a file documenting cell annotation is absent from their repository.We therefore reanalysed their data.Briefly, we mapped the raw reads to the hg38 reference genome.Seurat 21 (version 4.2.1) and Signac 22 (version 1.7.0) were used.We identified differential accessible regions of the clusters and signature genes for each cluster through overlapping analysis at promoters.We applied GO enrichment to the signature genes to identify overrepresented pathways for each cluster and used this information to annotate cell types from Li et al. 20 We then used the RunHarmony function in the Harmony package 23 to align our snATAC-Seq datasets with the public scATAC-Seq dataset, visualised cells in UMAP, projected cell type annotation from its respective data set and visually examined the positional relationship.

Tumour sample snATAC-Seq dataset extraction and annotation
After annotating the benign and tumour snATAC-Seq dataset, we extracted the tumour samples.For this extracted new object, Seurat 21 (version 4.2.1) and Signac 22 (version 1.7.0) were used for downstream analysis.Find-Clusters was used to identify clusters of cells by a shared nearest neighbours modularity optimisation based on SLM clustering algorithms 24 with a resolution of .8, and we identified 19 subclusters.Each cluster was generally dominated by one cell type after transferring annotation from the benign and tumour samples.Clusters were designated based on the predominant cell type and merged if they shared the same cell type.

ChromVAR analysis of snATAC-Seq dataset
For the snATAC-Seq dataset, RunChromVAR 39 function in Signac 22 (version 1.7.0) was used to calculate the transcription factor enrichment score.Motif matrices were obtained from the JASPAR2020 library.FindAllMarkers was used to identify differentially enriched transcription factor motifs for each subcluster against all other subclusters, 'wilcox' was set as the test used and logfc.thresholdwas set at .25.

GO enrichment analysis
Metascape 40 was used for GO enrichment analysis with the complete Homo sapiens proteome as background.We selected the top 10 hits for visualisation.

Survival analysis
The top 10 genes (by p-value) specific to a myeloid subcluster identified from differential accessibility analysis were used for survival analysis in the Survival Genie application. 41The Survival Genie Analysis Type selection criteria were 'Gene-Based', 'GENE SETS' and 'Upload my own'.Additional settings were select program: 'TCGA', select dataset: 'TCGA-PAAD, TCGA-CHOL, TCGA-COAD or TCGA-STAD', select tumour: 'Primary', select groups: 'Median', select survival: 'Overall Survival', select tumor infiltrating lymphocytes (TILS): LM22.The survival curves, hazard ratios and p-values are reported.

Single nuclear analysis of open chromatin shows unique cell populations in benign pancreas and PDAC
Benign pancreatic (n = 4) and PDAC (n = 4) samples were obtained following pancreaticoduodenectomy for eight total patients.Patients had pathology in the pancreatic head, which represents the majority of PDAC cases, 42 and had no neoadjuvant chemotherapy, as neoadjuvant therapy may modulate myeloid cell populations. 2,11,43Patients were of similar sex and age, and pathology was confirmed by histopathology of the resected specimen listed in Figures S1 and S2.Single nuclei suspensions were made directly from frozen tissues.Debris and ruptured nuclei were removed via fluorescence activated cell sorting (FACS), and we tested the transposase reaction for ATAC-Seq to identify a nucleosome banding pattern (Figure S3).After determining the method preserved chromatin integrity, we proceeded with chromatin accessibility in single nuclei.Combining benign and tumour samples into one UMAP, we showed similarity in the benign pancreatic samples but heterogeneity between PDAC projections (Figure 1A) and characteristic of PDAC tumours. 9,11,15hile there was some overlap in projections between nuclei from benign specimens and PDAC, the benign and tumour samples generally separated, with 18 694 nuclei identified for the benign specimens and 11 143 nuclei identified for the PDAC specimens (Figure 1B).We annotated eight cell clusters in the benign and tumour combined UMAP through an integrative analysis with the RNA part of snMultiome data for tumour and adjacent control tissues collected from an additional patient (Figures S4-S6), of which the cluster were annotated by expression of literature-supported markers (Figure S4).Marker genes of each snATAC-Seq cluster were predicted based on differential promoter openness (Figure S5B; see Table S2 for a full list of the predicted markers).For each snATAC-Seq cluster, we annotated each cell type by cross-referring its open chromatin markers to RNA markers for each cell type identified from snMultiome RNA data, calculated overrepresentations, and used that information to transfer the cell annotation at cluster level (Figure S5C).Representative marker genes from each of the snATAC-Seq clusters were then validated by gene expression in the snMultiome RNA data (Figure S5D).We further validated the annotation through inspecting the positional relationship of the clusters from the snMultiome ATAC when aligned with those from the snATAC-Seq dataset (Figure S6).Similarly, we validated the annotation with a published scATAC-Seq dataset derived from freshly acquired nuclei suspensions (Figure S7). 20he annotations for the snATAC-Seq data are summarised in Figure 1C.The majority of nuclei (53.0%) clustered into the acinar cell group (Figure 1D), aligning with literature showing acinar cells constitute approximately 50% of cells in healthy pancreas.Benign samples are predominated by acinar and ADM cells, while myeloid cells are specific to tumour samples (Figure 1E,F).We conducted a Fisher's exact test on the count of each cell type in both benign and tumour samples, yielding p-values consistently less than 5.3e-25.These findings show the unique cellular composition of the PDAC TME compared to benign pancreas identified by snATAC-Seq.

Open chromatin provides insight into the regulome and cellular function in benign and PDAC specimens
Chromatin accessibility can be used to infer transcription factor activity, gene promoter activity and therefore predicted molecular functions.Utilising chromVAR, 39 we identified differential binding of transcription factors across cell clusters (Figure 2A).In this analysis, acinar cells were enriched with T-cell factor/lymphoid enhancer factor 3 and 12 (TCF3 and TCF12) and Nuclear Receptor Subfamily 5 Group A Member 2 (NR5A2), transcription factors responsible for development and maintenance of the exocrine pancreas. 45In the ADM cluster, members of the TEA domain (TEAD) transcription factor family (TEAD1-4) have increased activity compared to the acinar cluster.TEAD transcription factors target SRY-Box Transcription Factor 9 (SOX9), 46 which drives ADM. 26 The ductal cell transcription factor activity was primarily driven by malignant ductal cells characterised by FOS/JUN activity, which comprise the activator protein 1 (AP-1) transcription factor 47 and have been implicated in PDAC development. 48ndocrine cell transcription factor activity was characterised by the regulatory factor X (RFX) transcription factor, important for development and maintenance of pancreatic beta cell function. 49Endothelial cell development is driven through ETS (E26 transformation specific) and ETV (ETS variant) transcription factors 50 also supported from the ChromVAR analysis.The fibroblast cluster showed activation of SMAD2/3 transcription factor activity, which is important for cancer-associated fibroblasts development and gemcitabine resistance in PDAC. 51The T cell cluster was also characterised by ETS transcription factors 1 and 2. ETS1 and ETS2 have differential functions in human T cells.ETS1 is high in resting T cells, while ETS2 mRNA and protein increase after T-cell activation. 52Transcription factor SPI1 is elevated in the myeloid subset as this is the primary transcription factor driving myeloid lineage cells and promotes monopoiesis over granulopoiesis and lymphoid cell development. 53Overall, ChromVAR analysis provides important transcription factor activity data consistent with cell type annotations, again confirming the validity of our cell type annotation for the snATAC-Seq data sets.
Next, we focused on genes with promoters preferentially opened in one cell type (Table S2).The promoter openings of those genes were first visualised using heatmap (Figure 2B).We next focused on promoter openings at genes known to be markers for each cell type (Figure 2C): regenerating family member 1 alpha (REG1A) is expressed by acinar cells. 11The ADM group was characterised by transient receptor potential vanilloid member 6 (TRPV6) promoter opening, implicated in promoting pancreatic cancer. 54Ductal cells in PDAC are characterised as keratin 19 (KRT19) positive. 11,15Endocrine cells have openings at the promoter for Paired box 6 (PAX6), important for development and maintenance of the endocrine pancreas in murine models. 27Promoter opening at the platelet and endothelial adhesion molecule 1 (PECAM1) promoter is characteristic of endothelial cells. 11Platelet-derived growth factor receptor beta (PDGFRB) is characteristic of PDAC fibroblasts. 15The myeloid population had opening at the Mannose receptor C-type 1 (MRC1) gene, characteristic of suppressive macrophages prevalent in PDAC. 5,12The T-cell population was characterised by promoter opening at CD69, an early activation marker for T cells expressed by T cells in PDAC. 55hen, we performed GO enrichment using the genes to identify unique molecular functions between the benign pancreas and PDAC.In the benign specimens, the GO enrichment was predominantly related to the nervous system, including axonogenesis, cell morphogenesis involved in neuron differentiation and axon development (Figure 2D), similar to previous analysis. 7These findings demonstrate the strong link between the nervous system and both the endocrine (alpha, beta and delta cells) and exocrine (acinar and ductal cells) pancreas, which are under control of the autonomic nervous system. 56In PDAC specimens, the hallmarks of the nervous system-pancreas crosstalk are replaced with functional terms for inflammation (positive regulation of cytokine-mediated signalling pathway), confirming previous scRNA-Seq analyses 15 and regulation of actin filament length, indicative of migratory or dividing cells (Figure 2E). 57Previous studies have also shown similar dysregulation of axonal guidance pathways in the progression to PDAC. 7e further applied GO enrichment to identify overrepresented molecular pathways in each cluster (Figure 2F; full list is provided in Table S4).Analysis of the acinar cell and endocrine clusters reveals the connection between the autonomic nervous system and pancreas for normal pancreatic development and function. 56The acinar cell cluster GO enrichment of regulation of neuron projection development contains the gene GFI1, which is known for its function in development of acinar cells. 58Furthermore, FGF9 and FGF10, genes recently shown to be involved in pancreatic development 59 and previously known to be important for neuronal development and maturation 60 play a role in the GO pathway developmental growth for the acinar cell cluster.In the endocrine cluster, the GO pathway neuron projection morphogenesis is characterised by Pax6, a gene responsible for development of the pancreas and other neuroectodermal structures, 61 ISL1, a gene important for beta cell function, 62 and ITGB1, a gene coding for an integrin critical for islet cell development in the pancreas 63 and for neuronal cell migration. 64DM GO enrichment includes cell-cell junctions and adhesions, representing the transition of acinar cells to a more ductal-like phenotype. 15In the ADM GO enrichment, PCDH1, a negative prognostic indicator in PDAC 65,66 is one of the genes for homophilic cell adhesion via plasma membrane adhesion molecules.Also, interestingly, the gene HEG1, which codes for a membrane protein that may be targeted in mesothelioma, 67 plays a role in the GO enrichment for cell-cell junction organisation.The ductal cells, which were primarily composed of tumour cells, have GO enrichment related to actin filaments. 57pecifically, the genes CAPG 68 and VASP, 69 which are related to metastasis in gastrointestinal cancers are represented in the regulation of actin polymerisation GO biological process identified as a function of the ductal cells.The endothelial cells function in blood vessel and vasculature development.In the endothelial cell cluster, the angiogenesis GO pathway was characterised by genes NOS3, important in VEGF-induced angiogenesis 70 and implicated as a therapeutic target in PDAC, 71 and FGF1, a known inducer of angiogenesis. 72The fibroblast cluster showed expected functions of extracellular matrix and structure organisation, 11 which were characterised by collagen genes (COL1A2, COL5A1, COL8A1, COL16A1).Interestingly, the fibroblast GO enrichment showed functions of regulating cell morphogenesis.Cancer-associated fibroblasts can activate epithelial to mesenchymal transition (EMT) in pancreatic cancer, 73 and the fibroblast GO biological pathway was also characterised by SPARC expression, which promotes tumour cell invasiveness in triple negative breast cancer when released by fibroblasts. 74As previously shown, 75 T cell activation involves regulation of GTPase mediated signal transduction.The gene ARHGAP15 was involved in the GO biological process for T cells, and previous studies have shown that different isoforms of ARHGAP15 may play different roles in T cell activation. 76For the myeloid cluster GO enrichment, activation of immune response was characterised by the inflammasomes (NLRP1, NRLP3 and AIM2), which can activate inflammation in myeloid cells 77 but may cause immunosuppressive T cells in PDAC. 78This function of myeloid activation of the immune system may be contrary to previous reports of myeloid cells in PDAC having immunosuppressive effects. 2,79Interestingly, the myeloid GO biological process, cellular response to cytokine stimulus, was characterised by genes MRC1 and CSF2RA, which are commonly associated with more immunosuppressive phenotypes in myeloid cells, 5,11,12,38 demonstrating that the myeloid cells responding to cytokine stimuli may be reprogrammed to an immunosuppressive phenotype.These findings suggest a spectrum of myeloid effector states in the PDAC TME, which should be further stratified.Taken together, these findings represent snATAC-Seq as a powerful tool for cellular heterogeneity delineation, regulome characterisation and molecular function inference from frozen pancreatic tissue.

High myeloid cell activity is characteristic of PDAC tumours but does not predict patient survival
Given the importance of myeloid cells in PDAC and inherent challenges with characterising these highly plastic cells, we next sought to determine if snATAC-Seq could give insight into myeloid cells.First, we determined that ChromVAR-inferred binding of transcription factors important for driving myeloid cell development (such as SPI1 and CEBPA) were more significantly more prevalent in tumour than benign pancreata (Figures 3A,B).These data further demonstrate the findings in Figure 1, in which only 24 nuclei representing myeloid cells (1.6% of the total myeloid population) were identified in the benign pancreata, consistent with previous analyses showing limited myeloid cells in normal pancreas. 12,20Given, the low number of myeloid nuclei identified in the benign pancreata, we further characterised the myeloid subpopulations in the PDAC TME alone.To do this, we re-clustered only the PDAC samples and annotated all the cell populations.We then focused on myeloid cells in the tumour samples (Figure 3C) and identified the top 10 genes with differential promoter openings in this cluster relative to other clusters (Figure 3D).Stratifying patients by high or low myeloid gene signatures using the Survival Genie application with the pancreatic adenocarcinoma (PAAD) dataset from The Cancer Genome Atlas (TCGA), we found no difference in overall survival between patient populations based on high or low expression of the tumour myeloid cell gene signature (Figure 3D).Since previous studies showed differences in disease progression or patient prognosis based on myeloid subsets in the tumour, 2,5,9,79 we next sought to stratify the myeloid cells.

Unique myeloid cell populations exist in PDAC tumour and have distinct functions
Observing no survival difference between groups with high and low myeloid signature gene expression, we further delineated the myeloid cluster into subclusters corresponding to several types of myeloid cells. 11We extracted the myeloid cell populations from the PDAC tumour samples (Figure 4A), performed re-clustering analysis and identified nine subclusters (Figure 4B).Most cells grouped into subclusters 0-4, which had contributions from all the assayed tumours (Figure 4B).The four smaller clusters of nuclei (subclusters 5-8) comprising less than 50 nuclei each and not represented in all tumours were excluded from further analysis (Figure 4B).Our snATAC-Seq analysis revealed each patient contributed differently to each myeloid subcluster (Figure 4C) consistent with the concept of tumour immune cell heterogeneity. 9,11Overall, analy-sis of myeloid cell populations from frozen PDAC tumours identified five major myeloid cell populations for further analysis type and function.
To further characterise the myeloid subpopulations, we determined their unique promoter openings, transcription factor activity and functional profiles.We first identified promoter regions preferentially open in each subcluster (Figure 5A).We noted no promoter upregulation in subcluster 0, suggesting the cells might be in a transitional state.In myeloid subcluster 1, we identified promoters for monocyte-related genes (ALOX5 and NLRC4).Arachidonate 5-lipoxygenase (ALOX5) is important for generating leukotrienes from arachidonic acid, 31 while NLR family CARD domain containing 4 (NLRC4) is a monocyte marker 32 important for activation of the innate immune response. 80Myeloid subcluster 2 is characterised by semaphorin 4A (SEMA4A), a marker for dendritic cell activation of T cells, 33,34 and semaphorin 7A (SEMA7A), which promotes dendritic cell migration. 35Myeloid subcluster 3 was identified as MDSCs based on S100 calcium binding protein A9 (S100A9), colony stimulating factor 2 receptor subunit alpha (CSF2RA) and interferon gamma receptor 2 (IFNGR2) expression. 11,38yeloid subcluster 4 was characterised by collagen expression (COL16A1, COL4A3, COL22A1, COL18A1) and Yes1-associated transcription regulator (YAP1), all TRM markers in murine PDAC. 5 Myeloid subcluster 4 also expresses keratin 7 and 8 (KRT7 and KRT8) aligning with epithelial-like myeloid cells identified by scRNA-Seq in PDAC.11 ChromVAR analysis of the transcription factor activity between the myeloid clusters was mostly distinct with some overlap as expected from highly plastic myeloid populations (Figures 5B,C).Myeloid subcluster 0 is characterised by FOXO transcription factor activity, which is important in maintenance of haematopoietic stem cells 81 and may be important in peritoneal TRMs.82 These cells may be a transitional state of myeloid cell as pancreatic TRMs can be replaced by stem cells upon damage to the pancreas.83 Myeloid subclusters 1 (recruited monocytes) and 2 (dendritic cells) were characterised by IRF1 F I G U R E 5 Characterisation of myeloid subclusters in pancreatic adenocarcinoma (PDAC) tumours.(A) Heatmap of the assay for transposase accessible chromatin (ATAC) signal at promoters showing cluster-specific chromatin accessibility within the myeloid subclusters. Te signal was averaged across all cells within each subcluster.Represented genes were shown on the left (p < .01).(B) Heatmap of top transcription factor activity, a promoter of anti-tumour immunity in conventional dendritic cells.84 IRF8 binding is also higher in myeloid subclusters 1 (recruited monocytes) and 2 (dendritic cells) and promotes differentiation of monocytes and dendritic cells.85,86 In myeloid subcluster 1 (recruited monocytes), the SPIC transcription factor, promoting anti-inflammatory function in monocytes 87 and angiogenic tumour-associated macrophages, 88 is characteristic.In contrast, FLI1 transcription factor in myeloid subcluster 2 (dendritic cells) is elevated, which drives development of phagocytic monocytes in mice. 89eloid subcluster 3 (MDSCs) has higher JunB and lower IRF8 activity.While JunB can drive macrophage activation, 90 it also drives MDSCs.91 IRF8 activity, which suppresses MDSC development, is lower in myeloid subcluster 3 (MDSCs).92 The transcription factor activity for myeloid subcluster 4 (tissue resident/epitheliallike macrophages) is similar to the activity for ADM, resembling the epithelial-like myeloid cells shown by scRNA-Seq in PDAC.11 Myeloid subcluster 4 also has elevated FOXO4 activity, whose activity has been shown to downregulate genes in peritoneal TRMs after feeding.82 Furthermore, NFIX, elevated in myeloid subcluster 4 (tissue resident/epithelial-like macrophages), can drive fibrosis in macrophages in murine models of muscular dystrophy.93 Taken together, these findings give insight into specific regulome of each myeloid subset.

Genotypic signatures of dendritic cells and tissue-resident macrophages may predict survival in PDAC patients
Having annotated the identity and functions of the myeloid subclusters, we next endeavored to determine if these myeloid subpopulations impacted survival.To assess a large patient cohort, we explored gene signatures for myeloid subpopulations in the TCGA-PAAD dataset.Using the Survival Genie 41 application, the signature genes for myeloid subcluster 1 (recruited monocytes) showed no survival difference, likely due to the mixed nature of this subset, where classical (pro-inflammatory) and nonclassical (anti-inflammatory) monocytes may have disparate effects on survival (Figure 6A).In contrast, high expression of the myeloid subcluster 2 (dendritic cell) signature in the tumour improved patient survival in the TCGA-PAAD cohort (Figure 6B), which has been similarly reported for circulating dendritic cells. 36We demonstrate no difference in survival based on the myeloid subcluster 3 (MDSC) gene signature (Figure 6C).While not statistically significant in the TCGA-PAAD cohort, there is a trend (p = .063)towards improved survival in patients with lower gene signature for myeloid subcluster 4 (tissue resident/epitheliallike macrophages) (Figure 6D).We then compared the significant gene signatures for myeloid subclusters 2 (dendritic cells) and 4 (tissue resident/epithelial-like macrophages) in other gastrointestinal malignancies, including colon adenocarcinoma (TCGA-COAD), cholangiocarcinoma (TCGA-CHOL) and stomach adenocarcinoma (TCGA-STAD) and saw no significant differences in survival (Figure S8).These findings support previous findings demonstrating that the prognostic significance of the myeloid cell infiltrate is unique to PDAC. 11ince we identified significant improvement in survival for the cluster 2 signature (dendritic cells) and a trend towards significant detriment to survival for the cluster 4 signature (tissue resident/epithelial-like macrophages) in the TCGA-PAAD dataset, we performed fluorescent immunohistochemistry to qualitatively identify costaining of signature genes that may be important for the  function of these cells, SEMA4A and TRIM29, respectively, with known myeloid markers.SEMA4A, which improves the ability of dendritic cells to activate T-cell responses 33,34 and was elevated in our gene signatures for subcluster 2, was co-stained with CD11c, a dendritic cell marker (Figure 6E).TRIM29 is implicated in the immunosuppressive functions of macrophages, 96 and we also identified TRIM29, which can function in the nucleoplasm, 97 expression in the nuclei of CD11b-expressing myeloid cells (Figure 6F).These data qualitatively validate the expression of SEMA4A and TRIM29 with known myeloid cell markers, CD11c and CD11b, respectively.Taken together, these data indicate that snATAC-Seq analysis can provide signature genes of interest for myeloid cells, such as infil-trating dendritic cells (myeloid subcluster 2) and tissue resident/epithelial-like macrophages (myeloid subcluster 4), which can be further explored in the preclinical and clinical settings.

DISCUSSION
In this study, we demonstrate snATAC-Seq can determine the cell compositions of frozen benign pancreas and PDAC samples.Assaying the accessible chromatin regions provides rich information on regulome landscape including transcription factor activities and gene promoter openness.These data can be utilised to determine cell type and infer molecular functionality.In PDAC specifically, previous studies have demonstrated that ATAC-Seq can be used to determine differentially accessible chromatin peaks in Epcam + epithelial cells sorted from surgically resected PDAC specimens.These peaks can then be used to predict patient disease-free survival. 98Another study utilising scATAC-Seq from freshly acquired benign pancreatic, normal adjacent tissue and PDAC specimens demonstrated that the KRAS mutation locks open an enhancer network in human PDAC that is similar to that found in murine PDAC models. 20Furthermore, another technique to look at open chromatin regions (single cell COOL-seq) has been used in PDAC specimens to demonstrate specific transcription factor binding motifs that may be important in pancreatic cancer progression. 7To add to this body of literature, we employed snATAC-Seq to explore the intricate regulome landscape of the diverse myeloid populations present in PDAC tumours.We demonstrated the composition differences among cell populations presented in benign pancreata and PDAC and validated our cell type annotations with snMultiomics and previously published scATAC-Seq, 20 adding rigour to our approach.snATAC-Seq was also used to compare the function of the benign pancreas, related to the pancreatic connection to the autonomic nervous system, 56 to the function of the cells in the PDAC tumour, which is related to infiltration of inflammatory cells similar to previous reports. 15Our GO enrichment determined the function of individual cell populations with support from the literature: the acinar and endocrine cells had links to the nervous system, 1,56 due to their functions related to neuron projection.The fibroblasts were shown to regulate cell morphogenesis such as previous studies showing fibroblasts drive EMT in PDAC tumour cells. 73We similarly identified functions of myeloid subsets, including immune response activation and regulation of lymphocyte and leukocyte activation.These findings were counter to previous studies showing myeloid subsets are primarily immunosuppressive in PDAC, 2,11,79 prompting further inquiry.
While the current study used four patient tumour samples to characterise myeloid subpopulations, the limited number of patients precluded associations with prognosis.Therefore, we used the determined gene signatures within a larger dataset to determine the influence of the ATAC defined unique populations in the myeloid compartment on survival.Unlike previous studies of infiltrating monocytes, we saw no correlation to overall survival based on expression of signature genes from monocyte 2 or MDSC 79 subclusters.Importantly, the previous studies measured monocytes 2 and MDSCs 79 in circulation, which may explain the differences in findings.The gene signature for dendritic cells predicted better overall sur-vival, similar to previous reports for circulating dendritic cells. 36Previous studies also show infiltrating dendritic cells promote a T cell response in PDAC, 99 and our gene signature showed expression of SEMA4A, which promotes T-cell activation. 33,34Furthermore, there was a trend towards worse overall survival in patients with a tissue resident/epithelial-like macrophages gene signature, supporting studies showing TRMs responsible for PDAC progression in murine models 5 and chemotherapy resistance. 43Similar to previous studies, 5,11 the tissue resident/epithelial-like signature had higher collagen expression.We also uniquely demonstrated TRIM29 expression in the tissue resident/epithelial-like signature, which can drive macrophage immunosuppressive features. 96The TRM/epithelial-like population was also like the epithelial-like macrophage population shown by scRNA-Seq of PDAC, 11 suggesting that TRMs and epithelial-like macrophages may arise from a similar cell population in PDAC tumours.The diverse populations of myeloid cells shown in this study are similar to scRNA-Seq studies showing unique populations within the myeloid cluster. 11linical trials focused on taming myeloid cells in PDAC have only been modestly successful.CCR2 antagonism therapy, which prohibits recruitment of circulating monocytes into the tumour, improved tumour control based on RECIST criteria when combined with standard of care FOLFIRINOX. 3 However, combining CCR2 antagonism with gemcitabine/nab-paclitaxel did not improve outcomes and was countered by risk of lung toxicity. 4While the CCR2/CCL2 axis may recruit myeloid cells that are hijacked by the tumour to become immunosuppressive, 2,3 these cells can promote an anti-tumour response when activated to a pro-inflammatory state. 95Furthermore, studies in murine models have shown TRMs, not monocytes, are the problematic population in PDAC, 5 but these cells have not been targeted clinically.Additional clinical trials have utilised CD40 agonism to promote antigen presenting cells to stimulate a T cell response, 100 again with limited efficacy.This may be in part due to dendritic cells also having an immunosuppressive function in tumours by expressing cytotoxic T lymphocyte-associated protein 4 (CTLA-4) 101 or programmed cell death protein 1 (PD-1), 102 preventing antigen presentation and T cell activation.Thus, while dendritic cells are critical for T cell activation, 99,37 overactivation can be detrimental and cause immunosuppression. 37,101,102Such observations point to the heterogenous nature of the myeloid cell infiltrate in PDAC, which we demonstrate here.This heterogeneity is driven by the interplay between the many factors comprising the PDAC TME, such as tumour cells, stromal cells, immune cells and soluble factors, which can differently impact myeloid populations.For example, CSF1-CSF1R interactions between tumour cells and macrophages, 12 respectively, are important for activation of macrophages.Importantly, colony stimulating factor 1 (CSF1) has been shown to be important for maintaining TRM populations in PDAC, which are pro-fibrotic and pro-tumour, 5 and may be responsible for therapeutic resistance to standard of care chemotherapy gemcitabine. 43PDAC tumour cells can also inhibit development of dendritic cells by release of granulocyte colony stimulating factor (G-CSF), which inhibits IRF8-mediated development of dendritic cells in the bone marrow. 36Furthermore, pancreatic stellate cells from the PDAC TME drive infiltrating myeloid cells towards a MDSC state via STAT3 signalling, which can in turn inhibit T-cell function. 103In contrast, the combination of C-C motif chemokine ligand 2 (CCL2) and interferon gamma (IFN-g) can drive an anti-tumour phenotype in PDAC myeloid cells, and this is particularly important in the presence of CD40 agonist treatment. 95hus, technologies such as snATAC-Seq can be useful for understanding the presence and function of myeloid cell subpopulations, which we have shown in the present manuscript.Such an understanding could lead to a more targeted approach towards implementation of myeloid therapies rather than a 'copy/paste' or 'one size fits all' approach, which has failed to translate immunotherapies that have been successful in other cancers to PDAC. 104Further studies utilising snATAC-Seq may uncover additional myeloid targets for future preclinical and clinical studies, emphasising a bedside to bench to bedside approach.
The current study has several potential limitations.While the current manuscript offers methodology for analysing myeloid populations in PDAC tumour specimens across multiple institutions with diverse patient demographics and treatment regimens, the current study only utilises four patients with benign pancreatic pathologies and four patients with PDAC.Given the limited number of PDAC patients, additional studies are warranted to better understand the myeloid populations, as this limited sample size may not be sufficient to conclusively characterise all potential myeloid subpopulations and their functions given tumour and patient heterogeneity.However, the present method provides a means for assaying frozen samples, which may broaden the available patient cohort by allowing for assessment of banked samples.

CONCLUSION
Despite a small sample size used in our analysis, this manuscript demonstrates the strength of utilising frozen tumour samples with snATAC-Seq.Frozen samples can be used to study patients with known outcomes and treatment conditions to tease apart the intricacies of the PDAC TME, particularly the myeloid subset, which can be challenging to study given the plasticity of the cells.Additional studies are warranted to understand how neoadjuvant therapies can affect myeloid subpopulations and how those myeloid subpopulations may in turn affect therapy response.With the implementation of such studies, patient-specific myeloid targets could be identified to better utilise myeloid therapies in particular patient populations or to implement novel myeloid targets.

F
I G U R E 1 Cellular heterogeneity and unique cell population to pancreatic adenocarcinoma (PDAC).(A) Uniform manifold approximation and projection (UMAP) visualisation of benign and PDAC single nuclear assay for transposase accessible chromatin (ATAC) sequencing (snATAC-Seq) dataset.(B) UMAP visualisation comparing benign and PDAC.(C) UMAP visualisation of the cell clusters and their annotations.(D) Cell type distribution among benign and PDAC snATAC-Seq datasets.(E) Cell type distribution among each dataset of benign and PDAC samples, highlighting myeloid (purple) as a unique cell population presented in PDAC samples.(F) Cell distribution comparing benign and tumour contribution to cell types.

F I G U R E 2
Regulome and molecular characterisation of cell populations identified in benign and pancreatic adenocarcinoma (PDAC) single nuclear assay for transposase accessible chromatin (ATAC) sequencing (snATAC-Seq) datasets.(A) Heatmap of top enriched transcription factor motifs activity (row) across cells organised into annotated clusters (column).(B) Heatmap of the ATAC signal at promoters of genes showing cluster-specific chromatin accessibility.The signal was averaged across all cells within each cluster.Represented genes were shown on the left.(C) Genome browser tracks for ATAC signal across representative genes of each cell type (shown at the top) with differentially accessible regions (including the promoters) highlighted by orange strips.(D) Gene Ontology (GO) enrichment analysis results for genes with promoters differentially accessible in benign samples as compared to tumour samples.(E) GO enrichment analysis results of genes with promoters differentially accessible in PDAC tumour samples as compared to benign.(F) GO enrichment analysis results for genes with promoter specifically accessible in designated cell types (shown at the top) as compared to all others.GO enrichment p-value is indicated by bar color with scale showing the -log p-value to the right.nGenes represents the number of genes under this GO term.

F I G U R E 3
Myeloid populations are higher in pancreatic adenocarcinoma (PDAC) tumour than benign pancreas.(A) Violin plot for transcription factor activity corresponding to SPI1 motif for benign and PDAC tumours from ChromVAR analysis on single nuclear assay for transposase accessible chromatin (ATAC) sequencing (snATAC-Seq) data.The mean comparison p-value was calculated with wilcox.testmethod.(B) Violin plot for CEBPA motif.(C) Re-clustering and uniform manifold approximation and projection (UMAP) visualisation of cell type among PDAC tumours datasets.(D) Kaplan-Meier survival curve for 'The Cancer Genome Atlas (TCGA)-pancreatic adenocarcinoma (PAAD)' patients sorted by overall expression of myeloid signature genes identified from snATAC-Seq analysis.

F I G U R E 4
Myeloid subclusters in pancreatic adenocarcinoma (PDAC) tumours.(A) Uniform manifold approximation and projection (UMAP) visualisation of the PDAC tumour dataset only for the myeloid populations.(B) UMAP visualisation of the identified subclusters.(C) Distribution of cells from each PDAC tumour to the subclusters identified in (B).
enriched transcription factor motifs activity (row) across cells organised into subclusters from myeloid cells (column).(C) The violin plot showing ChromVAR-inferred transcription factor activity score for selected motifs across subclusters.(D-G) Gene Ontology (GO) enrichment analysis results for genes with promoters preferentially open in myeloid subcluster of 1 (D), 2 (E), 3 (F) and 4 (G).GO enrichment p-value is indicated by bar color with scale showing the -log p-value to the right.nGenes represents the number of genes under this GO term.